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Abstract 

Genetic drift and non-random mating seldom influence species with large breeding populations and high dispersal 
potential, characterized by unstructured gene pool and panmixia at a scale lower than the minimum dispersal range of 
individuals. In the present study, a set of nine microsatellite markers was developed and used to investigate the spatio- 
temporal genetic patterns of the holoplanktonic jellyfish Pelagia noctiluca (Scyphozoa) in the Southern Tyrrhenian Sea. 
Homozygote excess was detected at eight loci, and individuals exhibited intra-population relatedness higher than expected 
by chance in at least three samples. This result was supported by the presence of siblings in at least 5 out 8 samples, 4 of 
which contained full-sib in addition to half-sib dyads. Having tested and ruled out alternative explanations as null alleles, our 
results suggest the influence of reproductive and behavioural features in shaping the genetic structure of P. noctiluca, as 
outcomes of population genetics analyses pointed out. Indeed, the genetic differentiation among populations was globally 
small but highlighted: a) a spatial genetic patchiness uncorrelated with distance between sampling locations, and b) a 
significant genetic heterogeneity between samples collected in the same locations in different years. Therefore, despite its 
extreme dispersal potential, P. noctiluca does not maintain a single homogenous population, but rather these jellyfish 
appear to have intra-bloom localized recruitment and/or individual cohesiveness, whereby siblings more likely swarm 
together as a single group and remain close after spawning events. These findings provide the first evidence of family 
structures and consequent genetic patchiness in a species with highly dispersive potential throughout its whole life cycle, 
contributing to understanding the patterns of dispersal and connectivity in marine environments. 
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Introduction 

Short life cycles, conferring the ability to mobilize huge 
populations that rapidly appear and disappear, allow gelatinous 
zooplankton to swiftly colonize available ecological spaces [1]. In 
particular, scyphozoans often form massive outbreaks that 
influence both ecosystem functioning and human activities 
(reviewed in [2,3]). 

Many medusozoans have complex life cycles with the succession 
of three stages: a short-living larva (planula), that metamorphoses 
into a benthic asexually reproducing polyp stage, giving rise to 
planktonic adult medusae through lateral budding (Hydrozoa), 
strobilation (Scyphozoa) or complete metamorphosis (Cubozoa). 
The presence of benthic polyps may limit the distribution of 
pelagic jellyfish in coastal areas with available hard substrates for 
larval settiement and polyp formation [1]. On the other hand, the 
modularity of a post-larval stage (the polyp) leads to exponential 
growth of jellyfish populations via polyembryony, i.e. the multiple 



production of adult (medusa) stages from a single fertilized egg [4] . 
However, some jellyfish species such as the mauve stinger Pelagia 
noctiluca are holoplanktonic (i.e. do not have a benthic polyp stage), 
and have a wide distribution in both inshore and offshore waters 
[5]. 

P. noctiluca has a global distribution [6] and inhabits preferen- 
tially warm and temperate waters [7-9]. In the Mediterranean 
Sea, massive blooms of P. noctiluca have been reported since the 
seventeenth century [10], insomuch that the species is almost 
considered a plague for human activities in coastal waters and has 
attracted special interest and concern since the late seventies [1 1]. 
Aggregations of hundreds thousands, or millions of individuals are 
not rare for this species whose population abundances show large 
fluctuations. Indeed, sudden demographical outbreaks lasting one 
or more years are normally followed by abrupt collapses, 
apparently with a periodicity of about 12 years [10,12]. The 
mechanisms leading to such oscillations are not completely clear. 
Goy et al. [10] suggested that P. noctiluca could take advantage 
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from occasional warm and dry weather conditions during late 
spring, going toward a demographic explosion in the following 
summer. This explanation is reasonable, as warm conditions could 
support the survival of a bigger number of the ephyrae produced 
by external fertilization of the eggs during spring-summer period 
[13]. Nevertheless, climatic conditions do not probably represent 
the only factor in play and further efforts are required to clarify the 
eco-physiological optima for P. noctiluca, identify biological 
characters and habitat changes apt to boost its outbreaks, 
distinguish critical thresholds of change, and quantify spatial and 
temporal levels of probabilities of outbreaks. By identification of 
reproductive units and measurements of gene flows, population 
genetics can help to fill some of current gaps of knowledge and 
advance our understanding of the ecological processes behind 
jellyfish bloom events. 

Previous population genetics and phylogeographic studies 
focusing on P. noctiluca at a large scale in the Mediterranean Sea 
and Eastern Atlantic Ocean indicated that high level of gene flow 
allows for a great connectivity across very large areas, maintaining 
a substantial panmixia [14,15]. Nevertheless, the peculiar popu- 
lation dynamics of this species in the Mediterranean Sea and in the 
neighbouring European Atlantic waters suggest the existence of 
complex underlying mechanisms. Indeed, several factors can play 
a substantial role in structuring natural populations, resulting 
sometimes in genetic patterns more complex than expected [16- 
18]. For instance, due to their high census population size many 
marine animal populations are generally considered not signifi- 
cantly influenced by processes as genetic drift or non-random 
mating. Moreover, when high population size is accompanied by 
great dispersal potential, unstructured panmictic populations are 
expected at a scale lower than the minimum distance dispersal of 
individuals. Nevertheless, a growing number of studies focusing on 
species characterised by high dispersal potential, as benthic marine 
invertebrates with a pelagic larval stage, showed low but significant 
levels of genetic differentiation among samples taken at distances 
far below the expected potential dispersal range [19-26]. Most of 
these studies highlighted also a co-occurrence between genetic 
heterogeneity on a small scale and temporal instability of genetic 
differentiation among populations, reporting swift fluctuations in 
time in the form of spatial structure changes across generations, or 
changes in allelic frequencies at a given sampling point. Such 
complex spatio-temporal genetic patterns led Johnson and Black 
[18] to coin the term "chaotic genetic patchiness" (CGP), later 
paraphrased by David et al. [27] as "fluctuating genetic mosaics". 
Different hypotheses have been formulated to explain CGP, but 
the most widely accepted is that factors linked to the reproductive 
strategy (e.g. high fecundity and high mortality in early life stages) 
can lead to a big variance in reproductive success and, 
consequendy, to a reduction of the effective population size 
[28,29]. These findings are changing a traditional paradigm, 
introducing the novel concept that many marine species, even if 
characterised by high census population sizes, can be exposed to 
processes usually considered effective only on small populations 
[30]. Hedgecock [28] compared the reproduction of such species 
to a sweepstakes lottery, characterised by an unbalanced 
distribution of the jackpot: "a small number of big winners grab 
all the prizes, leaving many losers empty-handed". Projecting this 
concept to a biological context, it means that stochastic factors can 
lead to an unbalanced genetic composition of the recruits due to a 
small number of successful progenitors. Consequence of this kind 
of processes could be a high level of relatedness between 
individuals coming from the same area and therefore inbreeding 
rates higher than expected. Nevertheless, even if many evidences 
support the existence of CGP processes, few studies have clearly 



demonstrated a connection between kin aggregation and "fluctu- 
ating genetic mosaics" patterns [23,31]. The aim of the present 
work is to use a set of newly developed microsatellite markers to 
study the population genetics of P. noctiluca at the small spatial scale 
of the Southern Tyrrhenian Sea and in a temporal range of 3 years 
(2010-2012), in order to elucidate the processes behind the genetic 
composition of each "Woom unit". Indeed, as already mentioned 
above, this species shows complex population dynamics that may 
suggest the driving action of unknown factors. Moreover, even if P. 
noctiluca is a holopelagic species, from a genetic point of view it may 
behave as a typical CGP species: it is characterised by high 
dispersal ability and the individuals of each bloom unit presumably 
spend together most part of their life, similarly to those belonging 
to a spatially stable benthic population. This particular habit 
suggests the possibility that related individuals, born at the same 
time and in the same area, could later spawn close to each other, 
producing inbred offspring and a bloom unit structure with high 
level of relatedness. Considering this hypothesis, the expected 
genetic patchiness could be further exacerbated by the high 
mortality rate due to stochastic factors, as unfavourable environ- 
mental conditions (i.e. sea storms), or by the high variance in 
reproduction success due to the external fertilization strategy 
[32,33]. The findings of the present work can play a pivotal role 
for any future investigations to explain the population genetic 
structure of the species at a broader geographical scale. 

Materials and Methods 

Ethics Statement 

Pelagia noctiluca is not an endangered species and no special 
permits were needed for sampling. All the animals were released 
without serious damages after sampling. Cnidarians are also 
renowned for their high regeneration potential. 

Study sites and samples collection 

A total of 259 individuals, belonging to the species Pelagia 
noctiluca (Forsskal 1775), were collected from 4 different locations in 
the Southern Tyrrhenian Sea, and from one site in the Northern 
Adriatic Sea (see Figure 1 and Table 1, for details on number of 
individuals analysed for each location, year of collection and 
coordinates). 

The collection sites in the Southern Tyrrhenian Sea are rather 
close to each other (the linear distance between the extreme point 
of the area is 215 km) and connected by a superficial eastward 
current [34,35] likely driving the passive transportation of 
planktonic organisms across the whole area. 

In order to address the species' genetic patterns on both spatial 
and temporal scale, the sampling was replicated in three 
consecutive years (2010, 2011, 2012) for Ustica (UST) and in 
two consecutive years (2011, 2012) for Messina (MES). No 
temporal replicates were available for the Northern Adriatic Sea 
(NAD06), Ischia (ISC 10), and Lipari (LIP11) locations, respec- 
tively sampled only in 2006, 2010 and 2011. The samples were 
collected by hand-net during bloom events choosing exclusively 
adults of the same size class, ranging from 7 to 10 cm of bell 
diameter. For each individual, a small piece of bell margin or oral 
arm was excised, preserved in 95% ethanol and stored at — 20°C 
until DNA extraction. 

Microsatellite markers development 

Total genomic DNA was extracted from the ethanol preserved 
tissues, following a CTAB-phenol-chloroform based protocol 
[36,37]. 
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Northern Adriatic Sea 
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Figure 1. Map of the study locations. Squares indicate locations for which temporal replicates are available; circles indicate location without 

temporal replicates. 

doi:1 0.1 371 /journal.pone.0099647.g001 



Microsatellite sequences were isolated by Ecogenics GmbH 
Switzerland (www.ecogenics.ch) using the high-throughput geno- 
mic sequencing approach described by Abdelkrim et al. [38], Two 
|J.g of genomic DNA from 1 2 P. noctiluca individuals belonging to 
the only two populations already available in 2010, Ustica (UST) 
and Ischia (ISC), were analysed on a Roche 454 GS-FLX platform 
(Roche, Switzerland) using 1/8 run and the GS-FLX titanium 
reagents. The total 53,066 reads had an average length of 277 base 
pairs. Of these, 843 contained a microsatellite insert with a tetra- 
or a trinucleotide of at least 6 repeat units or a dinucleotide of at 
least 10 repeat units. Suitable primer design was possible for 122 
reads, 1 2 of which were tested in agarose by Ecogenics GmbH 
Switzerland. In order to complement the set of loci, 1 0 additional 
reads were selected from the microsatellite library provided by 
Ecogenics using the software MSATCOMMANDER 1.0.8-beta 
[39] and primer pairs were designed with Primer3 v.0.4.0 (http:// 
bioinfo.ut.ee/primer3-0. 4.0/) [40,41]. In total, 22 primer pairs 
were tested on 20 P. noctiluca individuals of a single population 
(UST 10), performing PCRs as follows. The PCR volume of 10 ul 
contained approximately 50 ng of genomic DNA, 1 x Mastermix 
(RBC Taq Polymerase Kit), 0.5 u,M of a fluorescently labelled (6- 
FAM) universal Ml 3 primer (5 ' C AC G AC GTTGTAAAAC - 
GAC3') and the species-specific reverse primer, 0.15 u.M of 
species-specific forward primer with a 5' Ml 3 tail [42], 0.2 mM 



Table 1. Pelagia noctiluca sampling information. 



each dNTP and 1 unit taq. Amplifications were performed 
separately for each locus in a Eppendorf Mastercycler Gradient 
using the following thermal cycling profile: 94°C for 10 min, 10 
cycles of 94°C for 30 sec, 55°C for 45 sec, 72°C for 45 sec, 25 
cycles of 94°C for 30 sec, 53°C for 45 sec, 72°C for 45 sec and 
72°C for 10 min. Products were capillary electrophoresed on an 
ABI 3730XLs by Macrogen (Korea - http:/ /www.macrogen.com) 
using the internal size standard LIZ 500 (Applied Biosistems). The 
chromatograms were analysed using GeneMarker v. 2.2.0 (Soft- 
Genetics) and a final panel of nine loci was obtained after 
excluding loci with low signal, unclear peaks, excessive stuttering 
and the non-polymorphic ones. 

Markers genotyping 

Two multiplex panels were created with the help of the software 
Multiplex Manager 1.0 [43] setting 5 as maximum number of loci 
per reaction, 6 as primers complementarity threshold and 40 bp as 
minimum distance between loci of the same dye colour. Forward 
primers were labelled with different fluorescent dyes (6-FAM, 
VIC, NED - Applied Biosystems) and PCRs were carried out by 
multiplexing between 4 and 5 loci (multiplex content in Table 2) in 
1 0 (J.1 total volume containing: 5 ju.1 of Qiagen Multiplex PCR 
Master Mix, 2 |LLl of Q solution (Qiagen), 0.2 |J.M of each primer, 
and 50 ng of genomic DNA. PCR profile followed manufacturer's 



Sampling Site 


Sampling Year 


Abbreviation 


N 


Co-ordinates 


Northern Adriatic Sea 


2006 


NAD06 


52 


44°23'N, 14 44'E 


Ustica Island 


2010 


UST10 


43 


38°41'N, 1310'E 


Ustica Island 


2011 


UST11 


36 


38°41'N, 13 10'E 


Ustica Island 


2012 


UST12 


24 


38°41'N, 1310'E 


Ischia Island 


2010 


ISC10 


13 


40°44'N, 13'56'E 


Lipari Island 


2011 


LIP11 


53 


38°28'N, 14°57'E 


Messina Strait 


2011 


MES11 


14 


38 D 13'N, 15=38'E 


Messina Strait 


2012 


MES12 


24 


38°13'N, 1S=3S'E 


The table displays sampling sites and co 


ection years, the population 


sample acronym, the number of individuals and the 


sampling coordinates. 
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instructions and annealing temperature was set at 57°C for each 
multiplex reaction. The products of each multiplex PCR were 
electrophoretically separated on an ABI 3730XLs by Macrogen 
(Korea- http://www.macrogen.com) using the internal size 
standard LIZ 500 (ABI). Re-extraction and repeated multiplex 
amplifications were performed on individuals with failed PCRs. In 
addition, to assess loci amplification and scoring repeatability, 
10% of total sampled individuals were randomly re-amplified and 
alleles sized [44]. Allele sizes were assigned using GeneMarker v. 
2.2.0 (SoftGenetics) and allele scoring was repeated independendy 
by two authors /operators and then compared. Binning was 
automated with the software Flexibin ver. 2 [45] . All input files for 
further analysis were produced with the software Create [46] . 



m m Ln m m m 



<f- r- r- rM 



Markers characteristics and within population genetic 
variation 

Levels of polymorphism were calculated for each locus in each 
population. Number of alleles (Na) was obtained using the 
software FSTAT v. 2.9.3.2 [47], while observed and expected 
heterozygosity (H 0 ; H E ) were calculated by GENETIX v. 4.05.2 
[48]. 

Microsatellite loci are informative when they are independent 
one another and do not produce redundant information, therefore 
linkage disequilibrium was tested for each pair of loci within each 
population using Fisher's exact tests by Genepop web version 
[49,50] (http://genepop.curtin.edu.au/). The corrected signifi- 
cance threshold for multiple tests was set using the Bonferroni 
procedure [51,52]. 

As mentioned, the dataset was carefully checked against 
inconsistencies due to human mistake, genotyping errors and 
amplification failure. However, in order to evaluate and minimise 
the presence of artefacts, the final dataset was analysed with 
Microchecker v. 2.2.3 [53] and null allele frequencies estimated 
with the correction algorithm of van Oosterhout et al. [53]. Null 
alleles are a common feature of the microsatellites markers and are 
often accountable for an increased observed homozygosity. The 
software adjusts the number of homozygote genotypes in each size 
class to reflect the estimated frequency of null alleles and the 
"real" number of homozygotes. Therefore, a new dataset of each 
locus was obtained by considering the adjusted genotypes and, 
whenever possible, it was used simultaneously with the original 
one. The effect of the possible presence of null alleles was 
evaluated comparing the results obtained using the two datasets. 

In order to verify the effect of the correction for null alleles on 
HWE, two parallel analyses (with and without correction for null 
alleles) were performed with Arlequin ver. 3.5.1.3 [54], and 
respective p-values of HWE were calculated implementing an 
exact test with 1,000 steps in Markov chain and 10,000 
dememorization steps. Moreover, the same strategy was used to 
calculate the Weir and Cockerham [55] inbreeding indices (Fis) 
using FSTAT [47]. Significance levels for the global F IS values 
were calculated by performing 1 000 randomizations of genotypes 
among samples. 

High Fjs values and deviation from HWE expectations may be 
interpreted as an evidence of null alleles presence. Nevertheless, 
the Hardy-Weinberg disequilibrium (HWD) can be also due to 
phenomena as Wahlund effect [56] or non-random mating, 
especially when the HWD is due to heterozygotes deficiency and is 
associated with highly positive inbreeding values. The comparative 
approach with and without null alleles correction employed in this 
study aims to exclude null alleles pervasiveness. 
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Relatedness patterns 

The degree of relatedness between individuals of the same 
"bloom unit" was firstly investigated using the software ML-Relate 
[57]. Remarkably, this software allows testing for the presence of 
null alleles, as indicated by a deficiency of heterozygotes relative to 
Hardy-Weinberg expectations [58], and eventually accommodat- 
ing them in the subsequent analyses of relatedness. ML-Relate was 
therefore used employing a maximum likelihood approach to 
calculate pairwise Wright's [59] coefficients of relatedness (r) and 
the specific patterns of relationship (R) between individuals, 
classified as unrelated (UR), half-sib (HS) or full-sib (FS) and 
ranked depending on their likelihood values [ML(R)] . In order to 
verify if the presumed presence of null alleles increased the 
estimated relatedness between individuals within and among 
population samples pairs (inter- and intra-population r), two 
parallel analyses were performed with and without correction for 
null alleles. Relatedness values were obtained for all individual 
pairs (with and without accounting for null alleles) for the whole 
dataset of 259 specimens and used to calculate an average 
relatedness coefficient for each population pair. A two tailed t-test 
was then applied to verify if r-values obtained without correcting 
for null alleles were significantly different than those calculated 
taking null alleles into account. In addition, to verify if, by chance, 
intra-population r-values could have been equal or higher than 
observed, a Monte Carlo Markov Chain (MCMC) randomization 
of r coefficients was also applied calculating new inter- and intra- 
population r, after 1000 permutation of individuals across 
populations using PopTools [60] . Finally, a ranking of the most 
likely relationship was calculated for each individual pair within 
each population sample, and the frequency of pairs resulting to be 
FS and HS was estimated. Frequencies were calculated for both 
corrected and uncorrected datasets and compared by means of a 
two tailed t-test. 

Since every method of parental analysis and sibship reconstruc- 
tion includes a certain level of uncertainty ([61] and references 
therein), additional analyses were performed with the software 
Colony v. 2.0.4.1 [62]. Indeed, comparing the results of 
independent analyses performed with different programs is a 
strategy that can highlight possible failings in case of discordant 
results or ensure more robustness to agreeing outcomes. The 
statistical significance of differences between the two software's 
outcomes was assessed using a two-tailed t-test. 

The maximum full-likelihood method implemented in Colony 
v. 2.0.4.1, was used to partition dyads as HS or FS and to infer 
family structures within the "bloom units". Since Colony 
associates a probability value to each result, only those with a 
probability higher or equal to 95% were taken into account. 
Separate runs were performed for each sample, using the 
inbreeding model, as suggested for dioecious species when the 
inbreeding level is high (cf. Colony user guide), and setting 
medium length of the run and high likelihood precision. Due to 
the lack of prior information about sibship size, the value 0 was set 
up for known male and female genotypes as well as for known 
paternal and maternal sibship. Since the program Colony can 
accommodate null alleles and other stochastic genotyping errors in 
the analysis, once again two runs were performed for each 
sampling site, taking or not the presence of null alleles into 
account. In the runs "without nulls", the genotyping error was set 
to the default level (0.005), whereas the runs "with nulls" were 
performed allowing locus specific genotyping error levels, obtained 
adding the null alleles frequencies calculated by Microchecker v. 
2.2.3 for each locus in each population to the default genotyping 
error rate. HS and FS dyads frequencies were calculated for both 



corrected and uncorrected datasets and then compared by means 
of a two tailed t-test. 

Family structures were inferred reconstructing extended sibship 
networks, namely clusters of individuals connected by a chain of 
HS or FS intermediate individuals. In practice, if the individual A 
shares a parent with the individual B, and the individual B shares a 
parent with the individual C, in that case A and C are linked 
through B and are members of the same extended sibship network, 
even if they do not have a common parent. Since the extended 
sibship networks can also include FS, they were represented as FS 
families (groups of individuals sharing both the parents) nested into 
sibship networks. For each sibship network the number of FS 
families (and the respective number of family members) was 
reported, when any. 

Population structure 

To examine patterns of genetic variation among the studied 
population samples, pairwise F ST values and corresponding p- 
values were calculated with Arlequin ver. 3.5.1.3. Prior to this 
analysis, in order to estimate the extent of bias possibly introduced 
by the presence of null alleles, pairwise Fg-p were also calculated 
using the software FreeNA [63] with and without null alleles 
correction. This software, indeed, estimates null allele frequencies 
for each locus and population, following the Expectation 
Maximization (EM) algorithm of Dempster et al. [64] and then 
implements the so-called ENA correction to provide accurate 
estimation of Fgx [65] in presence of null alleles. 95% confidence 
intervals for the FgT values were obtained using 50,000 bootstrap 
iterations. F S x estimates obtained with and without applying the 
ENA algorithm were compared by means of a two tailed t-test. 

Results 

Markers characteristics and within population genetic 
variation 

Of 22 initially selected loci, 13 were not polymorphic or showed 
very low signal, not clear peaks or excessive stuttering and were 
discarded, while 9 markers were suitable for the present study and 
resulted to be polymorphic in all population samples (Table 2). 
The total number of alleles ranged from 5 to 18 (Mean = 8, 
Standard Deviation SD = 3.7) while observed heterozygosity 
ranged from 0.14 to 0.85 (Mean = 0.48, Standard Deviation 
SD = 0.16) (Table 3). 

Microchecker did not identify scoring errors associated with 
stuttering, but suggested the presence of null alleles by analysing 
strong deviations from HWE (Table 3). Indeed, all loci showed 
heterozygote deficit in at least 3 out of 8 populations, except for 
the locus Pelnoc_16756 having a moderate homozygote excess 
only in one population (LIP1 1, P = 0.030). Loci potentially affected 
by null alleles over all populations are showed in Table SI. 

Two out of 36 pair-wise locus comparisons revealed significant 
linkage disequilibrium, after standard Bonferroni adjustment 
[51,52], between loci Pelnoc_46263 vs. Pelnoc_40199 and 
Pelnoc_46263 vs. Pelnoc_7445 at UST10 only. This result might 
be due to lack of recombination, non-random mating in inbred 
populations or admixture of genetically distinct populations (i.e. 
Wahlund effect [56]), all factors able to cause loci to appear 
statistically linked [66] . 

All loci but one (Pelnoc_ 16756; P-value = 0.290) did not globally 
result in HWE. HWD was due to heterozygote deficiency for all 
cases, suggesting two alternative explanations: the presence of 
technical artefacts, such as non-amplifying alleles (null alleles), or 
the influence of biological factors [67]. However, given the strong 
deviation from HWE observed in almost the whole set of loci, null 



PLOS ONE | www.plosone.org 



5 



June 2014 | Volume 9 | Issue 6 | e99647 



Inbreeding, Relatedness and Chaotic Patchiness in Pelagia noctiluca 



o 
o 
o 
d 
V 
o. 



I 



v 

a. 



o ■- 

V 9 
a. o 



V 
o. 



x 



a 3 



yo co co m 



v£> O r\| m CO O t- 

vo <— o co 10 
oo m co m co m io 



0 
X 



z 



o 
Q 
< 

z 



0 
X 



ro vo CTi m 



x 



o ,_ 



V V V 



o 

V 



vO 00 00 



CTi ^- i- O 

m io io 



co oo m co m co 



< 
Z 



\D OO CO CO 



m r\i m ro 



^ ' — ' — vO O"* CO < — 



0 
I 



z 



m m io yo 



o 

V 



o o 

V V 



x 



V 

0. 



V 
a. 



m fo cri 

in rv i— m 
co ro co ro 



vo oo in co 



< — in m ' — 



0 
X 



ci a* cn co 
cn m in m 



0 
X 



< 
z 



z 



io in co 



X 



v ^ 



o 
o 
o 



CO ro ai a* 



CO CO vo 



oo oo m co 



CO i£> vO 



a 
o 
a. 



o 
X 



z 



ro yo m yo yo 



in vo m m o 
M m io o i- 



oo o> vo in 

in m in * 

3 i- * 

O O VO IV 

* * i- o 



a 
o 



o 
I 



z 



VO 0^ CO t- 



in \o m m 

in in if) o 

V0 * IN O 

O 51 10 T 



IM 
ST 
O 



o> in 3 

■-is * 

o 10 N 

* *- o 



aj 5 co 
S 1 S 

1 s§ 

c 1/1 aj 
„ J c 
ra o O 

y, U Q. 

2 ? ^ 
°> 2 5 

1/1 a> o 



a ai ai at 



PL0S ONE | www.plosone.org 



6 



June 2014 | Volume 9 | Issue 6 | e99647 



Inbreeding, Relatedness and Chaotic Patchiness in Pelagia noctiluca 



alleles specifically affecting each single locus seem unlikely to be 
the only factor responsible for the observed pattern of disequilib- 
rium. In addition, if null alleles only were involved, the HWD 
should disappear after correction using Microchecker (see 
methods). Indeed, although for some loci in some populations 
the extent of HWD for the corrected dataset was lower (HWE p- 
values are shown in Table S2), globally the HWE was not 
recovered. In fact, only for 1 4 out of 4 1 test showing HWD (locus 
by locus estimation) the equilibrium was re-established after 
applying the correction for null alleles. 

Inbreeding coefficient estimation confirmed the significant 
excess of homozygotes, with Fjs values generally higher than zero 
(Table 4). Similarly to what observed for the HWE comparisons, 
inbreeding estimates calculated using the dataset corrected for null 
alleles were still positive (average 0.087), although lower than those 
calculated with the original, uncorrected, dataset (average 0.301). 

The observed heterozygosity lower than expected at almost all 
loci, the presence of sporadic linkage disequilibrium in one single 
population, and the persistence of HWD and positive F IS values 
after correction for null alleles represent altogether strong 
complementary indications that true biological factors, rather 
than technical artefacts only, are responsible for the observed 
pattern of disequilibrium. 

Relatedness patterns 

Average relatedness between individuals within and among 
Southern Tyrrhenian population samples was investigated using 
ML-relate with and without correction for null alleles. The two- 
tailed t-test applied to compare r-values obtained with and without 
correction resulted to be significant (P<0.0001) due to higher r- 
values obtained taking into account the presence of null alleles 
(Relatedness Monte Carlo Simulations are shown in Tables S3 
and S4). However, a higher than expected relatedness among 
individuals of several samples was detected that was robust against 
null allele correction (Table 5). Indeed, Monte Carlo simulation 
tests indicate that 4 out of 8 intra-population pairwise comparisons 
resulted to have a significantly higher r-value than expected by 
chance when not taking into account null alleles presence (Table 5). 
Correcting for null alleles resulted in a small difference and 3 out 
of 8 significant comparisons remained significant after the 
correction. Considering conservatively the results obtained with 
null alleles correction, our approach suggested that at least 



UST10, MESH and MES12 have higher within-population 
relatedness degree than expected by chance. 

When using ML-Relate to estimate the parentage relationships 
among individuals (Table 6A), all samples, except ISC 10, MES1 1 
and UST12, were shown to contain some related individuals. The 
highest frequencies of HS and FS pairs were present in LIP 1 1 
(HS = 0.03544; FS = 0.00452) and in NAD06 (HS = 0.03619; 
FS = 0.00150). Moreover, the amount of HS and FS pairs 
calculated with and without accounting for null alleles was not 
significantly different (p-values>0.30), clearly indicating that this 
result cannot be explained by non-amplifying alleles. 

The results of the analyses of parentage relationships performed 
by Colony show agreement with those obtained using ML-Relate 
by disclosing the existence of related individuals in many of the 
investigated populations. Separate approaches are used by the two 
software: Colony jointly considers the likelihood of larger patterns 
of relationship, whereas ML-Relate independently determines the 
relationship of each pair of progeny [68]. Therefore, some 
differences were also found. In particular, HS dyads were found 
by Colony in all samples except ISC 10 and MES11, while FS 
dyads were found in UST10, LIP 11 MES12 and NAD06 
(Table 6B). Unlike the ML-Relate results, the highest HS dyads 
frequencies are present in UST10 (0.04858), while the highest FS 
pairs frequencies are shown by MES12 (0.00724). Also in this case, 
the comparison between relationships calculated with and without 
accounting for null alleles was not significant (P>0.18), whereas 
the two tailed t-test between the ML-Relate and Colony outcomes 
showed a significant difference only between HS frequencies 
(P = 0.028) due to a higher average of HS frequencies found by 
Colony (Mean = 0.031, Standard Deviation SD = 0.019) respect to 
ML-Relate (Mean = 0.0 13, SD = 0.015). 

As regards family reconstruction with Colony, results were 
consistent among replicates. Since the results of the analyses with 
and without accounting for null alleles were not significandy 
different, only the results of the runs "without nulls" are shown. 
The more structured families (containing both HS and FS 
individuals) were found in UST10, LIP11, MES12 and NAD06, 
with respectively 3 FS families nested within 7 sibship networks, 4 
FS families nested within 9 sibship networks, 2 FS families nested 
within 1 2 sibship networks and 5 FS families nested within 1 1 
sibship networks. All FS families were composed by two 
individuals. A graphic representation of the family structures 



Table 4. F| S values calculated with and without null alleles correction. 





FIS Estimates 


Locus 


Uncorrected 


95% C.I. 


Corrected 


95% C.I. 


Pelnoc 40622 


0.317 




-0.008 




Pelnoc_39456 


0.428 




0.091 




Pelnoc 46263 


0.318 




0.112 




Pelnoc_44003 


0.310 




0.122 




Pelnoc_44210 


0.380 




0.040 




Pelnoc 40428 


0.201 




0.100 




Pelnoc_40199 


0.414 




0.149 




Pelnoc_16756 


0.069 




n.e. 




Pelnoc_07445 


0.150 




0.051 




TOT 


0.301 


0.214-0.371 


0.087 


0.058-0.113 



Locus by locus and over all loci F !s values were calculated for the original and the corrected dataset after Microchecker v. 2.2.3 analysis, 
n.e.: Not evaluated due to HWE and no changes in the corrected dataset. 
doi:1 0.1 371 /joumal.pone.0099647.t004 
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Table 6. Frequencies of Half Sib and Fu 


I Sib dyads 


estimated by A) ML-Relate and B) Colony. 






A 


ruuuiaiiuiij 


Frequencies 












HS NNA 




HS NA 


FS NNA 


FS NA 


NAD06 


0.03619 




0.02488 


0.00150 


0,00301 


UST10 


0.01357 




0.01432 


0.00226 


0,00226 


UST11 


0.01131 




0.00980 


0.00000 


0,00075 


UST12 


0.00000 




0.00075 


0.00075 


0,00075 


ISC10 


0.00000 




0.00000 


0.00000 


0.00000 


LIP11 


0.03544 




0.03544 


0.00452 


0,0015 


MES11 


0.00000 




0.00000 


0.00000 


0.00000 


MES12 


0.00377 




0.00301 


0.00075 


0,00075 


p-values 


0.3250 






0.8436 




B 


Populations 


Frequencies 












HS_NNA 




HS_NA 


FS_NNA 


FS_NA 


NAD06 


0.04444 




0.06570 


0.00580 


0.00097 


UST10 


0.04858 




0.07152 


0.00539 


0,00404 


UST11 


0.03968 




0.06507 


0.00000 


0.00000 


UST12 


0.02766 




0.06719 


0.00000 


0.00000 


ISC10 


0.00000 




0.00000 


0.00000 


0.00000 


LIP11 


0.03338 




0.04789 


0.00290 


0,00072 


MES11 


0.00000 




0.01282 


0.00000 


0.00000 


MES12 


0.04710 




0.47101 


0.00724 


0,00724 


p-values 


0.2096 






0.1839 




HS: half-sib dyads; FS: full-sib dyads; NNA: no null alleles accounted for; NA: null alleles accounted for. 

In each table, the last row reports the p-values of the t test comparing each category of estimated frequencies (e.g. FS "with nulls' 
doi:1 0.1 371 /journal.pone.0099647.t006 


vs. FS "without nulls"). 



found in these 4 populations is shown in Figure 2. Relating to the 
other investigated populations, Colony identified in UST11 a 
single sibship network composed by 28 individuals, while other 7 
individuals resulted unrelated two sibship networks, respectively of 
10 and 5 individuals, were identified in UST12, paralleled by 8 
unrelated individuals. None of the latter two populations showed 
FS families, whereas the populations ISC 10 and MESH did not 
show any kind of family structure. 

Population structure 

Population pairwise Fgx values were initially calculated with 
FreeNA to verify if accounting for null alleles could bias this 
genetic distance estimation (Table S5). F ST values calculated with 
and without applying the ENA correction resulted to be 
significantly different due to an increase of F ST distances when 
accounting for null alleles. The occurrence of null alleles produces 
an overestimation of the Fst values in case of significant 
population differentiation, which is dampened by corrections for 
nulls [63]; therefore, only uncorrected F ST values were considered, 
believing they still represent the best estimations for a conservative 
analysis of population differentiation. Accordingly, the global and 
the pairwise F ST values calculated with Arlequin ver. 3.5.1.3 from 
the original dataset were taken into account for further consid- 
eration. The overall Fst value was small (0.01714, 95% CI 
0.00501-0.03949) but highly significant (P<0.0001). Inspection of 
pairwise F ST values (Table 7) indicated that the overall differen- 



tiation found is mainly attributable to the sample collected at 
Messina in 2012, which was significantly different in 5 out of the 7 
test performed. Surprisingly, this differentiation gready exceeds 
that of the sample from Northern Adriatic Sea (NAD06). Despite of 
more than a thousand kilometres of geographical separation from 
the Tyrrhenian samples and a different year of collection (2006 
against 2010-2012), they resulted to be significantly differentiated 
only in 2 out of 7 test, whereas all the other comparisons were not 
significant (P^0. 06934). Interestingly, the distinctiveness of the 
MES12 sample does not seem to reflect stable geographic 
differentiation, given that MES12 was different from the sample 
MES 1 1 , collected one year before at the same location (MES 1 1 - 
MES12 F ST = 0.05338, P<0.0001), and homogeneous with the 
sample UST11, collected in 2011 at —230 Km of distance 
(UST1 1-MES12 F ST = 0.00184, P = 0.90137). Moreover, tempo- 
ral variation was identified also at the Ustica site (UST 1 0-UST 1 2 
F ST = 0.04007, P<0.0001). 

Discussion 

Pelagia noctiluca blooming populations in the Southern Tyr- 
rhenian Sea exhibited significant deviation from HWE due to 
large excess of homozygotes for 8 out of 9 microsatellite loci, 
leading to high inbreeding coefficients (Fjg). Moreover, higher 
relatedness than expected by chance, between individuals within 
and among population samples was detected and supported by the 
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1(2) 




1(2) 



1 




1 






1(2) 


1(2) 


K2) 



Figure 2. Family structure in four Pelagia noctiluca Mediterranean populations. The first line of each dendrogram shows the extended 
sibship networks produced by Colony v. 2.0.4.1. For each population, in this line the squares indicate the number of individuals linked by kinship 
(directly or indirectly). For example the number "34"indicates that 34 individuals are linked together by kinship. Separate squares refer to individuals 
with no connections. The second line shows the Full Sib (FS) families nested in the sibship network. Namely, since each network includes all the 
individuals linked by a certain level of kinship (individuals sharing one or two parents), when any, in the second line the FS individuals (sharing both 
parents) were reported, specifying the number of families and the respective number of family members (in brackets). For example, the numbers 
"4(2)" indicate the presence in the upper sibship network of 4 FS families composed by 2 individuals each. 
doi:1 0.1 371 /journal.pone.0099647.g002 



presence of half-sib (HS) individuals in at least 5 out of 8 samples, 4 
of which (at least) contained full-sib (FS) individuals also. The 
genetic differentiation among samples was globally small, but it 
highlighted a spatial and temporal genetic patchiness probably 
reflecting the influence of reproductive processes, as also suggested 
by the outcomes of the relatedness analyses. 

In most studies of population genetics, large deviations from 
HWE are expected to be strongly linked to null alleles, and for this 
reason several markers are frequently discarded. Although the 
occurrence of some null alleles in microsatellite markers cannot be 
completely ruled out, however HWD can also have biological 
explanations, especially if the deviations are chronic at multiple 
loci [69,70]. In the present study, all loci but one strongly deviated 
from HWE expectations. Therefore, even considering the 
presence of null alleles, it is unlikely the only factor accountable 
for the HWD and the positive F IS values, as supported by the 
missing recovery of HWE (27 out of 41 locus by locus evaluations 
remained in HWD after null alleles correction, Table S2) and by 
the still positive F IS values after the null alleles correction (Table 4). 
Moreover, if null alleles number of blanks was observed in the loci 
Pelnoc_44210 and Pelnoc_39456 (respectively 9 and 8 blanks over 
259 individuals) while all the other loci showed a number of blanks 
ranging from 0 (Pelnoc_46263) to 5 (Pelnoc_07445) with an 
average number of 2 blanks over 259 individuals. According to the 
estimated frequencies of null alleles, for example the locus 



Pelnoc_46263 should show at least five blanks instead of the zero 
observed, as well as the locus Pelnoc_40199 showed only three 
instead of the expected nine blanks; on average 4 null homozy- 
gotes should be observed for each locus (data not shown). Basing 
on this line of reasoning and considering also the sporadic linkage 
disequilibrium observed in one single population (UST10) as a 
possible symptom of non-random mating in inbred populations 
[67], null alleles can be reasonably excluded as the only cause of 
the strong HWE. Conversely, available evidence conveys links 
between the observed HWE and biological factors, as supported 
by the high intra-population relatedness degree (Table 5) and by 
the occurrence of related individuals in several sampled popula- 
tions (Table 6), that were found both considering or not the 
presence of null alleles. 

P. noctiluca is a species with external fertilization that reproduces 
during spring and summer [71,72]. Broadcast spawning in the 
water column would theoretically ensure homogeneous distribu- 
tion of gametes and random fertilization, but reproductive and 
behavioural features may drive patterns of kin aggregation. This 
seems to be the case for P. noctiluca. First, each mature female 
jellyfish spawns oocytes in a sticky mucus ribbon, holding eggs 
together for several minutes before its dissolution [72]. This 
peculiarity may favour fertilization of the whole set of oocytes by 
sperms released by a single or a few male mates, producing a large 
amount of full sibs. Second, fusion of gametes produced by related 
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individuals and the resultant formation of inbred offspring may be 
favoured by aggregative swimming behaviour of jellyfish. Indeed, 
P. noctiluca en-masse proliferations can be characterised by high 
densities active aggregations routinely exceeding 100 individuals 
per cubic meter, driven as a whole by surface winds and marine 
currents [73,9,11]. Related individuals born from the same 
parental group at the same time have a reasonable probability 
to remain together in the native bloom unit also during subsequent 
spawning events [74], increasing the probability to produce inbred 
olfspring.Indeed, P. noctiluca jellyfish display an active swimming 
behaviour during daily vertical migration, in response to a 
corresponding circadian migration pattern of their zooplankton 
prey [75-77]. Canepa et al. [78] also suggested the occurrence of a 
seasonal migratory circuit along vertical corridors in the proximity 
of marine canyons, acting as circulation pumps of water bodies 
favouring local jellyfish aggregations. Physical oceanographic 
models suggest that water circulation is characterised by both 
downwelling and upwelling events, influencing nutrient exchange, 
biological productivity and eventually the composition of shallow 
and deep-sea biota [79-81]. Enhanced upwelling near canyon 
areas provides increased nutrient exchange that boosts phyto- 
plankton and, hence, zooplankton abundance. Available evidence 
suggests indeed that submarine canyons have important effects on 
coastal marine ecosystems, including food webs [79], by acting as 
additional drivers of environmental and biological discontinuities 
of the coastal habitats. Pelagia noctiluca aggregations are known as 
strongly influenced by marine currents and favoured by physical 
discontinuities as fronts and pycnoclines [73,78,82]. The up- and 
downwelling currents driven by marine canyons may therefore 
represent a driving force for dense aggregations, for both physical 
(jellyfish pooled together by water movements) and trophic (local 
increasing of preys availability) reasons [78]. 

Interestingly, an atypical swimming behaviour of P. noctiluca 
jellyfish has been also reported throughout spring and summer 
months [78], with frequent formation of jellyfish couples which 
may presumably boost fertilisation rates as well as full-sibs 
generation. Living in swarms can provide strong advantages to 
jellyfish not only facilitating conspecific gametes to fuse, but also 
allowing highly synchronised reproduction among conspecifics to 
enhance fertilisation success. Furthermore, although many marine 
reproductive cycles appear to be on lunar, circadian or circatidal 
rhythms, waterborne chemical cues are crucial for fine-tuning 
spawning synchrony [83,84]. 

Comparably, at least in some fish species, individuals can 
remain together from birth to settlement, despite relatively long 
planktonic durations. Sensorial and behavioural mechanisms 
enable fish larvae to remain in close proximity of each other 
throughout their planktonic dispersal and achieve genetically 
homogeneous recruitment [85]. Finally, individual dispersal in 
several marine taxa characterised by a pelagic larval phase, 
including sponges, echinoderms, molluscs, crustaceans, corals and 
fishes, may be influenced by oceanographic conditions [86-93] . In 
P. noctiluca the maintenance of kin-related jellyfish aggregation 
along marine currents may be similarly influenced by small-scale 
hydrodynamic and oceanographic patterns conducive to limited 
individual mixing despite high dispersal potential. Hydrogeo- 
graphic features such as eddies, gyres or upwelling fronts could 
restrict dispersal of groups of medusae in confined areas [82], 
preventing an extensive mixing with individuals belonging to other 
aggregations and allowing the fusion between gametes produced 
by related individuals. Consistently with this hypothesis, Lee et al. 
[94] recently suggested oceanographic barriers to dispersal causing 
genetic differentiation among some geographically near 
(—200 Km) populations of jellyfish PJiizostoma luteum in the Irish 
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Sea. Even if the bentho-pelagic life cycle of this species may favour 
the retention of medusae in coastal areas, a substantial genetic 
homogeneity should be expected at low spatial scale, unless to 
consider specific water circulation patterns as hindrance to gene 
flow. 

Additional evidences of biological factors influencing the P. 
noctiluca genetic structure come from population differentiation 
analyses. According to classical genetic theory, a marine species 
with high pelagic dispersal potential is expected to have no clear 
genetic structure at least below the minimum dispersal distance of 
individuals [95-98] and previous phylogeographic studies on P. 
noctiluca confirmed this expectation. Indeed, Miller et al. [15], in a 
phylogeographic study performed using mitochondrial cyto- 
chrome oxidase subunit I (COI) and two nuclear internal 
transcribed spacers (ITS1 and ITS2) genes, found high and 
statistically significant genetic differences (at two out of three 
markers: COI: 0> ST = 0.72, P<0.001; ITS2: 0 ST = 0.023, P< 
0.001) only between Southern and Northern Atlantic samples, 
geographically very far from each other (—10.000 km), whereas no 
genetic difference was found at a smaller geographic scale. The 
authors interpreted this result as a suggestion of historical rather 
than contemporary gene flow. Conversely, Stopar et al. [14] 
provided support for present-day panmixia among P. noctiluca 
populations from the Mediterranean Sea and the North-western 
Atlantic Ocean. These authors could not detect significant genetic 
differentiation relative to both COI and IST1/ITS2 markers, 
apart from differences between samples from Northern and 
Southern Adriatic Sea (COI: <1> ST = 0.095, P<0.()1). Such a short- 
scale pattern of genetic differentiation is not coherent with an 
isolation by distance model, and has been interpreted as the 
outcome of basin-scale hydrodynamic processes reducing the 
mixing of individuals born in different areas of the Adriatic Sea. 

Similarly, significant genetic differences were found in this study 
among populations just a few km far apart (i.e. Lipari-Messina 
— 70 Km; Lipari-Ustica = 150 Km; Ustica-Messina «230 Km). 
Conversely, 5 out of 7 pairwise comparisons between the 
populations from the Southern Tyrrhenian Sea and the one from 
the Adriatic Sea (NAD06, separated by a distance of more than 
1000 Km from the closest Tyrrhenian sample) did not result 
statistically significant (P^0. 06934). Moreover, pairwise compar- 
isons among consecutive samples collected in the same locations 
(Messina, Ustica) in different years showed significant genetic 
differences, highlighting the lack of temporal stability of the 
genetic composition of the P. noctiluca blooms (Table 7). Altogether 
these findings suggest the existence of biotic and abiotic 
mechanisms influencing the genetic pool of P. noctiluca to generate 
a temporally dynamic mosaic of small-scale genetically differen- 
tiated patches rather than a homogeneous mixing of the 
population, or a geographic set of populations isolated by distance. 

Genetic heterogeneity on a small geographic scale, especially 
when temporal genetic differentiation is stronger than spatial 
differentiation, is generally attributed to temporal changes in the 
genetic composition of recruits [19,21,27]. We suggest that the 
genetic structure observed in P. noctiluca can be explained as fine- 
scale genetic patchiness, which may be generated through 
processes driving localised temporal variation of numbers and 
genotypes of recruits [19]. In such a case, stochastic factors dealing 
with reproduction processes may influence the proportion of 
individuals contributing to the next generation, leading to 
temporal variance in allelic frequencies of the recruits [28]. Under 
such hypothesis, also called the "hypothesis of sweepstakes 
reproductive success" (SRS) [28], in species with high fecundity 
and high mortality rate at early stages, many individuals fail to 
contribute to recruitment. Several factors as local oceanographic 



conditions (such as occurrence of canyons and upwelling areas), 
short life-time of gametes [99] and temporary spatial constraint of 
individuals can act on fertilization success and formation of 
recruits pool, generating an "instantaneous drift effect" [27]. Due 
to the variable parental contribution to recruits pool, the genetic 
composition of recruits can change generation by generation, 
leading to spatio-temporal genetic patchiness [100]. From the 
evolutionary point of view, SRS has important consequences 
because, due to the stochastic nature of the process involved, 
divergence is not accumulating but renewed each generation 
[19,98,101] and can be counteracted in the long term by dispersal 
and gene flow [26,102]. On the other hand, potentially allowing 
the replacement of the entire population by a small fraction of 
individuals, SRS provides the power for rapid evolutionary change 
and for population resilience [103,104]. The SRS hypothesis also 
provides testable predictions [105] such as a reduction of effective 
population size, measurable effects of genetic drift though time, 
reduction of allelic diversity and increase of relatedness within 
cohorts. 

To date, many studies provided evidence of a widespread 
occurrence of chaotic genetic patchiness (CGP) in several 
benthic marine species [19-22,24-26]. However, only very few 
studies, on fishes [21,23,106] and barnacles [31], demonstrated 
a direct connection among intra population relatedness, family 
structure and CGP. By considering the high inbreeding and 
relatedness found at P. noctiluca population level, together with 
the presence of half- and full-siblings in several samples, the 
present study may contribute to clarify the role of family 
structure in CGP formation. The most probable scenario to 
explain the temporally unstable genetic patchiness of P. noctiluca 
populations is the co-occurrence of large variation in the 
reproductive success of individuals and genetic drift. Indeed, 
the high intra-population inbreeding level, highlighted by 
homozygote excess and positive Fi S values, combined with the 
presence of related individuals in several samples could be a 
symptom of variance in reproductive success generated by 
environmental, reproductive and behavioural factors. Moreover, 
the high mortality rate of ephyrae [107-109] could generate an 
instantaneous genetic drift able to enhance the formation of 
genetically unbalanced assemblage of recruits through random 
suppression of some allelic variants. 

Chaotic genetic variability among populations, explained by 
alternate periods of rarity and abundance and related genetic drift, 
is in accordance also with the "Flush and Crash" speciation model 
[110] as a driving microevolutionary force leading to genetic 
diversity also at small spatial and temporal scale [111]. P. noctiluca 
populations in the Southern Tyrrhenian Sea are seemingly 
characterised by a spatio-temporal CGP and random genetic drift 
may represent a leading micro-evolutionary force shaping the 
genetic structure of this species. Even if several studies have 
already shown similar patterns in species with high larval dispersal 
ability [32,112-116], this study provides the first evidence of 
family structures and consequent genetic patchiness in a highly 
dispersive holopelagic species. 
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Table S3 Relatedness Monte Carlo Simulation, null 
alleles accounted for (adjusted allele frequencies by ML- 
Relate), 1000 iterations. 
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Table S4 Relatedness Monte Carlo Simulation, no null 
alleles accounted for, 1000 iterations. 
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Table S5 Pairwise FST values calculated using the 
program FreeNA. The table S5A shows the results of the 
analyses performed using the uncorrected dataset. The table S5B 
is referred to the Fst values obtained implementing the ENA 
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